Avalanches and Waves in the Abelian Sandpile Model 
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We numerically study avalanches in the two dimensional Abelian sandpile model in terms of a 
sequence of waves of toppling events. Priezzhev et al [PRL 76, 2093 (1996)] have recently proposed 
exact results for the critical exponents in this model based on the existence of a proposed scaling 
relation for the difference in sizes of subsequent waves, As = Sk — Sk+i, where the size of the previous 
wave Sk was considered to be almost always an upper bound for the size of the next wave Sfc+i. 
Here we show that the significant contribution to As comes from waves that violate the bound; the 
average (As(s*,)) is actually negative and diverges with the system size, contradicting the proposed 
solution. 
PACS numbers: 64.60.Lx 



The sandpile was the first model introduced by Bak, 
Tang, and Wiesenfeld to demonstrate the principle of 
self-organized criticality Q . Self-organized criticality de- 
scribes a general property of slowly-driven dissipative 
systems with many degrees of freedom to evolve toward a 
stationary state where activity takes place intermittently 
in terms of bursts spanning all scales up to the system 
size. The sandpile model has subsequently received a 
great deal of attention due in part to its potential for 
having a theoretical solution. Dhar showed that certain 
aspects of its behavior could be calculated exactly based 
on the Abelian symmetry of topplings. For instance, rig- 
orous results have been obtained for the total number of 
allowed configurations on the attractor of recurrent states 
H , for some height-height correlation functions , and 
later for the distribution of sizes of the last wave in an 
avalanche ||, among other quantities. Nevertheless, a 
formal solution for the all important critical exponents 
describing the distribution of avalanche sizes and dura- 
tions in the Abelian sandpile model (ASM) has remained 
an elusive goal. 

In a recent Letter || exact results were proposed for 
the distribution of avalanche sizes in the ASM, based on a 
decomposition of an avalanche into a sequence of "waves" 
of topplings. Specifically, the prediction was made that 
the asymptotic distribution of the number of topplings in 
an avalanche is P(S) ~ S~ T with r = 6/5, and the dis- 
tribution of the number of sites covered by an avalanche 
is P(a) ~ a~ Ta with r a — 5/4. The results of numerical 
simulations do not convincingly support these claims, al- 
though measuring the actual avalanche distribution ex- 
ponents in this model is notoriously difficult. 

Here we scrutinize the main assumption in the argu- 
ment leading to the predicted exact results in Ref. || 
for the ASM. Based on careful numerical simulations 
we show that the fundamental assumption that the next 
wave usually is contained by the previous wave fails dras- 
tically. In particular the exponent a as defined in [pj does 
not exist, and the difference in sizes of subsequent waves 



As is more often negative than positive. In fact, the 
negative contribution has a sufficiently fat tail that the 
average difference (As(sfc)) is negative and diverges with 
system size for all Sk < s co . The quantity s co is the cut- 
off in waves sizes due to the finite size system, s co ~ L 2 , 
where L is the linear extent of the system. In short, the 
physics is dominated by waves that are not contained 
within their predecessors. Thus the argument leading to 
the claimed exact results is incorrect. 

The ASM consists of a square lattice of size L with 
a discrete number of sand grains occupying each site. 
Initially, the lattice may be empty, and sand is dropped 
grain by grain at random sites. After each drop, all sites 
which exceed a critical threshold for stability, Zi > z c = 3, 
are "toppled" by distributing a single grain of sand to 
each of their four nearest neighbors or, for boundary 
sites, over the edge of the lattice. Toppling proceeds for 
a number of time steps until all sites are stable again, 
and a new grain is dropped at a random site. Dropping 
sand represents an external driving force on the system 
whose impact is dissipated in intermittent sequences of 
toppling events which are called avalanches. The num- 
ber of toppling events following the addition of a single 
grain is the size S of an avalanche. Starting from an 
empty lattice, avalanches are initially rare and only of 
short duration. But the system fills up with sand to the 
point that many sites are close to threshold. Then, the 
system reaches a stationary state in which for any one 
grain dropped, on average, one grain must leave the sys- 
tem through the open boundaries. The grains are trans- 
ported by avalanches which are now broadly distributed 
in both duration and extent over many orders of magni- 
tude, only limited by the system's size. Thus, the sys- 
tem has self-organized into a critical (SOC) state with a 
highly correlated response to the external driving. 

This model has a few other remarkable properties, 
most notable the fact that the order of toppling events 
during an avalanche is interchangeable ( "Abelian" ) with- 
out changing the final state of the system, which for in- 
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stance allows an exact enumeration of the critical state 
. Also, it was found that the domain spanned by a sin- 
gle avalanche is always compact, though with a fractal 
boundary 0. Manna || introduced a different sandpile 
model, without Abelian symmetry, where the toppling 
grains are stochastically distributed to nearest neighbors 
so that the distribution is symmetric only on average. 
For some time it was believed, based on real space renor- 
malization group arguments H and Manna's numerical 
simulation results, that the ASM was in the same uni- 
versality class as the Manna model. But for both models 
the values of the distribution exponents, obtained by var- 
ious extrapolations of the results from extensive numer- 
ical simulations, have barely converged to within 10% 
after many years of study. Recently, Ben-Hur and Bi- 
ham computed the geometrical scaling properties of 
avalanches instead of individual avalanche distribution 
exponents. Their results indicate that the ASM and 
the Manna model may belong to different universality 
classes. A survey of two-time autocorrelation functions 
in various SOC models Jll| also reveals differences be- 
tween both models: the ASM exhibits "aging" |L2| while 
the Manna model does not 

Effort has focused recently on understanding avalanche 
dynamics in the ASM by decomposing the avalanche into 
a sequence of more elementary events. Dhar and Manna 
Q introduced the notion of inverse avalanches which 
were shown to be equivalent to a direct representation 
of avalanches in terms of waves of toppling events [ fhH . 
Ivashkevich, Ktitarev, and Priezzhev defined waves as 
follows: if the site (i) to which a grain was added be- 
comes unstable, topple it once and then topple all other 
sites that become unstable, keeping the initial site (i) 
from toppling a second time. The set of sites that top- 
pled thus far are called "the first wave of topplings" since 
every site can only topple once. After the first wave is 
completed the site (i) is allowed to topple the second 
time, not permitting it to topple again until the "second 
wave of topplings" is finished. The process continues un- 
til the site (i) becomes stable and the avalanche stops. 

This elegant decomposition of avalanches into waves 
reveals many interesting features. First of all, the waves 
are individually compact, and each site that topples in a 
wave topples exactly once in that wave. As a result, the 
state of the system after a wave is exactly the same as the 
state before the wave except at sites on the single closed 
boundary of the wave which separates sites that toppled 
in that wave from sites that did not. Just inside the 
wave boundary a trough relative to the previous heights 
appears, whereas just outside the boundary a hill appears 
where sand was transport outside of the wave. Thus the 
sequence of topplings in the next wave will follow exactly 
the sequence of topplings in the previous wave until the 
wave first reaches the prior wave's boundary, at which 
point the sequence may differ. Priezzhev et al || argued 
that generally subsequent waves are spatially contained 



within the previous waves because of the trough at the 
boundary. Their analysis only considers waves which are 
contained within previous waves and they "neglect the 
overlapping of waves and deal only with the decrease of 
wave size" ||. Using spanning tree arguments together 
with this assumption, they find that the size difference 
between subsequent waves < As >= Sk — Sk+i is pos- 
itive, finite in the limit of large system size, and obeys 
a scaling relation As ~ s£. Key to the argument lead- 
ing to Eq. (12) in Ref. Q is the length of the boundary 
r of the previous wave Sk which is presumed to contain 
the subsequent wave Sk+i- Here we show that the main 
contribution to As comes from waves which escape the 
boundary of their preceding wave, and are bounded only 
by the system size. In fact, the average (As(sfc)) is nega- 
tive and diverges to infinity as the system size increases. 

We have simulated about 10 7 waves in L 2 -systems up 
to L = 1024. To ensure the accuracy of our numerical 
simulations, we have reproduced a variety of previously 
obtained exact results for the distribution of waves. For 
instance, our data yields the s -11 / 8 power-law that was 
derived by Dhar and Manna Q for the distribution of 
the very last wave in each avalanche. We also found the 
1/s behavior for the distribution of all waves . 
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FIG. 1. Plot of the distribution P(sk+i\sk) for the next 
wave to be of size Sk+i, given that the previous wave was of 
size Sk in a system of size L = 1024. Each graph contains 
data for 2 m < s k < 2 m+1 for m = 4, 5, . . . , 19 from bottom to 
the top. To avoid overlaps, the graphs are offset. It appears 
that for each value of Sk, P(sk + i\sk) initially falls like s^{ 4 
(as the dashed line on top), but at a scale linear in Sk crosses 
over to a s^f{ 4 fall-off (as the lower dashed line). 

To determine As we recorded the distribution of sub- 
sequent waves for a given size of the preceding wave 
P(sk+i \sk) as shown in Fig. 1. In Fig. 2 we show the data 
collapse where the horizontal axis is now x — Sk+i/sk- 
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There are two regimes separated by a turning point near 
one. There is clear evidence of a fat power law tail for 
x 3> 1, which will dominate the average As for each 
value of Sk- The data is sufficiently well represented by 
a scaling form 



P{s k+1 \s k ) 



s~ F 



Sfc 



where F(x — > 0) — > 1 and F(x 3> 1) ~ x r 
collapse indicates that (3 ~ 3/4 and r ~ 1/2. 
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The data 
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FIG. 3. Plot of < As >= Sfc — Sfc+i as a function of 
the previous wave Sk, obtained from P(sfc+i|sfc) (see Fig. I) 
via Eq. ([]) for system size L — 2 8 , 2 9 and 2 10 from top to 
bottom. Initially, < As > is negative and falling for increas- 
ing s. Closer to the cut-off s co the linearly rising term in 
Eq. ([]) dominates. In each case, the graph turns positive at 
Sfc « L 2 /I0. 
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FIG. 2. Scaling collapse for the data in Fig. f according 

to Eq. (§), giving F(x) ~ s 4 +1 P(sfc + i|sfc) as a function of 
x = Sfc+i/sfc. The tail for each graph falls approximately like 
x~2 (as the dashed line). 

From Eq. (Q), the computation of (As(sfe)) =< Sk — 
Sfe+i > leads immediately to 



(As(sfe)) = Sfc - Sfc 



(2) 



where s co is the cutoff in wave sizes from the finite sys- 
tem size; s co ~ L 2 . For Sfc <C s co this gives < As(sfe) >= 



1 or, using our values for /3 and r, 



< As(sfc) >= sfc 



Oi co *fc 



(3) 



where C is a positive number that depends on the details 
of the function F. It is important to note here that (As) 
is negative for all Sfc up to a scale L 2 , and that it di- 
verges with the system size. Our numerical measurement 
for (As) shown in Fig. 3 confirms the above analysis. 

If we (erroneously) exclude from the data set all waves 
that are larger than their predecessor, i. e. eliminating 
those waves that contribute to a negative As, then we 
reproduce (in Fig. 4) the plot given in Fig. I of Ref. M. 
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FIG. 4. Plot of < As > as a function of the previous wave 
s, but leaving out all data where As is negative. We chose 
L = 512 to compare with Fig. f in Ref. ||. 

Furthermore, a more detailed analysis of subsequent 
waves that are smaller, in terms of number of topplings 
s, than their predecessor reveals that in most cases the 
following (smaller) wave still, more often than not, ex- 
ceeds the confines of the (larger) previous wave at some 
points on the boundary. The smaller waves actually do 
escape the boundary of their predecessor. As the system 
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size grows, the the fraction of consecutive waves which 
violate the assumption of Ref. also grows. 



In summary, we have shown that the analysis in Ref. 
H is fundamentally flawed because it deals only with 
waves of decreasing size, whereas the dominant contri- 
bution to As comes from the "fat tails" in the distribu- 
tion P(sk+i\sk) where waves escape the boundary of their 
predecessor explicitly. Yet, our numerical data seems to 
indicate that a certain degree of regularity in the distribu- 
tion of consecutive waves exists, leading to what appears 
to be exact values for f3 = 3/4 and r — 1/2. It may be 
possible that the spanning tree arguments used in Ref. 
H can be generalized to include the dominant contri- 
bution from overlapping waves. The appearance of the 
apparently simple exponents /3 and r gives us some hope 
that an exact solution using the elegant decomposition 
of avalanches into waves can be discovered. 
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